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Using density functional theory molecular dynamics simulations, we predict shock Hugoniot curves 
of precompressed methane up to 75 000 K for initial densities ranging from 0.35 to 0.70 gem -3 . At 
4000 K, we observe the transformation into a metallic, polymeric state consisting of long hydrocarbon 

chains. These chains persist when the sample is quenched to 300 K, leading to an increase in shock 
compression. At 6000 K, the sample transforms into a plasma composed of many, short-lived chemical 
species. We conclude by discussing implications for the interiors of Uranus and Neptune and 
analyzing the possibility of creating a superionic state of methane in high pressure experiments. 



Methane is one of the most abundant chemical species 
in the universe, with the vast majority of our solar sys- 
tem's share being locked away in the interiors of the gi- 
ant planets [T] . While recent work has suggested that icy 
components of the cores of gas giants such as Jupiter and 
Saturn are likely to be dissolved into the surrounding hy- 
drogen [2j |3] , the ice giants Uranus and Neptune likely 
contain a significant amount of methane at pressures 
ranging from 20 to 80 GPa and 2000 to 8000 K. Extraso- 
lar planets in the ice giant mass regime are now known to 
be extremely common throughout the universe @]. This 
motivates our need for a greater understanding of the be- 
haviour of this important material in this specific range 
of pressures and temperatures. 

Dynamic shock experiments [3J E| combined with static 
precompression [7HS] have the potential to probe the be- 
havior of methane at pressures and temperatures in ice 
giant interiors but these experiments are challenging and 
their interpretation often not straightforward [TO] HT] , 
There exists a long tradition of first-principles simulation 
predicting the outcome of past IT2T4T4"] and future [T5] 
shock experiments and providing an interpretation on the 
atomistic level for the thermodynamic changes that were 
observed |16j . Here we use density functional molecu- 
lar dynamics (DFT-MD) simulations to characterize the 
properties of dense, hot methane and identify a molecu- 
lar, a polymeric, and a plasma state. From the computed 
equation of state, we predict the shock Hugoniot curves 
for different initial densities in order to guide future shock 
wave experiments with precompressed samples. 

The high pressure, high temperature behavior of 
methane has previously been probed experimentally and 
theoretically across a range of conditions, revealing a 
complex chemical behavior characterized by the conver- 
sion of methane into higher hydrocarbons and eventu- 
ally diamond as pressure is increased. Laser heated dia- 
mond anvil cell (LHDAC) studies by Benedetti et al. [17] 
showed the presence of polymeric hydrocarbons and dia- 
mond at pressures between 10 and 50 GPa and tempera- 



tures of about 2000 to 3000 K. More recently, a LHDAC 
experiment by Hirai et al. [18] found evidence of the for- 
mation of carbon-carbon bonds occurring above 1100 K 
and 10 GPa, including the formation of diamond above 
3000 K. Gas gun shock wave experiments by Nellis et 
al. [H] achieved pressures of 20 to 60 GPa and produced 
results that were interpreted as the formation of carbon 
nanoparticles. Radousky et al. [20] previously carried out 
shock experiments at lower pressures that did not reach 
the dissociation regime. 

Ancilotto et al. [21] performed DFT-MD simulations 
that followed the predicted isentropes of the intermedi- 
ate layers of Uranus and Neptune. At 100 GPa and 4000 
K methane was found to dissociate into a mixture of hy- 
drocarbons. At 300 GPa and 5000 K it became apparent 
that carbon-carbon bonds were favored with increasing 
pressure. Recent DFT-MD simulations with free ener- 
gies derived from the velocity autocorrelation function 
by Spanu et al. [22] predicted that higher hydrocarbons 
than methane are energetically favored at temperature- 
pressure conditions between 1000 and 2000 K and at 4 
GPa and above. Static zero-kelvin calculations by Gao et 
al. |23j predicted methane to decompose into ethane at 95 
GPa, butane at 158 GPa and diamond at 287 GPa. Gold- 
man et al. [23] computed shock Hugoniot curves for the 
molecular regime of methane up to 50 GPa. Better agree- 
ment with experiments was obtained by approximately 
including quantum corrections to the classical motion. 
The vibrational spectrum for each atomic species was de- 
rived from DFT-MD simulations and the internal energy 
was corrected by assuming each atom has six harmonic 
degrees of freedom. This approach is exact for a har- 
monic solid but its validity remains to be more carefully 
evaluated for molecular fluids. 

Laser and magnetically driven shock experiments [S] |S] 
now allow multi-megabar pressures to be obtained rou- 
tinely, however the lack of independent control of the 
pressure and temperature variables makes the direct 
study of planetary interior conditions difficult [5] |H] ■ For 



2 



any given starting material, the pressure and tempera- 
ture conditions of the shocked material are uniquely de- 
termined by the shock velocity. Conservation of mass, 
momentum and energy across the shock front lead to the 
Hugoniot relation [2"5] . 

H = (E-E ) + ^(P + P )(V-V )=0, (1) 

relating the final values of the internal energy, E, pres- 
sure, P, and volume, V, to the corresponding values 
Eq,Pq, and Vo for the initial state of the material. While 
shocks launched into materials at ambient conditions al- 
low extremely high pressures to be reached, a large part 
comes from the thermal pressure and the density is often 
only about 4 times higher than the initial value |15j . The 
resulting shock temperatures are significantly higher than 
those along the planetary isentrope. Static precompres- 
sion of the sample, however, shifts the Hugoniot curves to 
higher densities allowing higher pressures to be explored 
for given temperature. Precompressed shock experiments 
on methane are thus highly desirable as a means of ex- 
ploring the pressure-temperature conditions of interest to 
planetary science. 

In this article, we present theoretically computed 
methane shock Hugoniot curves corresponding to initial 
densities ranging from 0.35 to 0.70 gcm~ 3 and temper- 
atures up to 75 000 K. We use DFT-MD simulations 
to determine equations of state of methane at relevant 
temperature-pressure conditions, compute the equation 
of state, and analyze the atomistic and electronic struc- 
ture of the hot, dense fluid in order to make predictions 
for future shock measurements using methane. 

COMPUTATIONAL METHODS 

All DFT-MD simulations in this article were per- 
formed with the Vienna Ab Initio Simulation Package 
(VASP) [33]. Wavefunctions were expanded in a plane 
wave basis with an energy cutoff of 900 eV, with core 
electrons represented by pseudopotentials of the projec- 
tor augmented wave type [27]. Exchange-correlation ef- 
fects were approximated with the functional of Perdew, 
Burke, and Ernzerhof [55]. Each simulation cell consisted 
of 27 methane molecules in a cubic cell with temperature 
controlled by a Nose-Hoover thermostat. The Brillouin 
zone was sampled with T point only in order to invest the 
available computer time into a simulation cell with more 
atoms. Electronic excitations were taken into account 
using Fermi-Dirac smearing |30j . 

Simulations were undertaken at 79 different 
temperature-density conditions as shown in Fig. [T] 
chosen in order to encompass the Hugoniot curves for 
methane at initial densities ranging from 0.35—0.70 
gcm~ 3 . Temperatures ranged from 300 to 75 000 K 
and pressures from 3 to 1100 GPa. Each simulation 
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FIG. 1: Methane shock Hugioniot curves showing temper- 
ature as a function of density for different initial densities, 
po = {0.35,0.40,0.45,0.50,0.55,0.60,0.65,0.70}. The sym- 
bols show the density and temperature conditions of our DFT- 
MD simulations that either stayed molecular, or reached a 
polymeric or plasma state. 

lasted between 1 and 2 ps, and in some cases up to 8 ps. 
Almost all simulations used an MD time step of 0.2 fs. 
Initially, we performed a few simulations with a time 
step of 0.4 fs at lower temperatures but reduced the 
timestep as we explore higher temperature and densities. 
Our tests showed that the computed thermodynamic 
functions at lower temperatures agreed within error bars 
with the results obtained with the shorter time step. 

RESULTS 

The pressure and internal energy computed from DFT- 
MD simulations at various densities and temperatures 
are given in Tab. [T] in the online supplemental informa- 
tion. Any initial transient behavior was removed from 
the trajectories before the thermodynamic averages were 
computed. Hugoniot curves were calculated by linear in- 
terpolation of the Hugoniot function, H, in Eq. |T]) as 
a function of temperature and pressure. Figure [T] shows 
temperature as a function of density along the Hugoniot 
curves generated for each of the eight initial densities, po- 
Several features are immediately apparent. As might be 
expected, precompression allows higher densities to be 
reached at any given temperature. Each Hugoniot curve 
can be divided into three segments. Up to 3000 K, one 
finds the temperature to increase linearly with density. 
In our simulations, the methane molecules remained in- 
tact. At a temperature of approximately 4000-5000 K, a 
plateau is reached where density increases rapidly with- 
out a corresponding increase in temperature. Our simu- 
lations show the system entering into a polymeric regime 
where the methane molecules spontaneously dissociate to 
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form long hydrocarbon chains that dissociate and reform 
rapidly. We will demonstrate this regime to be metallic, 
which implies a high electrical conductivity and reflectiv- 
ity should be detected during a shock experiments. 
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FIG. 2: The pressure as function of the compression ratio, 
pi pa, along the shock Hugoniot curves from Fig.[T]for different 
initial densities, po- 

In the third regime, the temperature on the Hugoniot 
curve again increases rapidly as the system assumes a 
plasma state where many ionic species exist for a very 
short time. The slope of the curve significantly depends 
on the initial density. For a small degree of precompres- 
sions, density on the Hugoniot curve does not increase in 
the plasma regime because the sample has already been 
4-fold compressed. For a high degree of precompres- 
sion, density of the Hugoniot curve keeps increasing in 
the plasma regime because the theoretical high tempera- 
ture limit of 4-fold compression has not yet been reached. 
This is confirmed in Fig. [2j where we plotted the shock 
pressure as a function of compression ratio, p/po- The 
highest compression ratio of 4.0 is obtained for 10 000 K 
and the lowest initial density. In general, the compression 
ratio is controlled by the excitation of internal degrees of 
freedom that increase the compression and interaction 
effects that act to reduce it [T5]. The strength of the in- 
teractions increases with density, which explains why the 
compression ratio decreases when the sample is precom- 
pressed. 

We now focus on characterizing the polymeric state 
that gives rise to the predicted increase in shock com- 
pression at 4000 K. Fig. [3] shows a linear increase of pres- 
sure with temperature as long as the methane molecules 
remain intact. When the molecules dissociate at 4000 K, 
the pressure reaches a plateau at 40 GPa for 1.2 gem -3 
before it rises gradually above 5000 K when the sys- 
tem turns into a plasma. For a higher density of 
1.5 gem -3 , one observes a significant drop in pressure 
when molecules disassociate that resembles the behavior 
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FIG. 3: The upper panel shows the pressure-temperature re- 
lation of CH4 DFT-MD simulations that have been heated 
(filled symbols) and of samples that have been cooled (open 
symbols) from a polymeric state at 4000 K. The lower panel 
shows the corresponding electronic gap bands and its vari- 
ance during the MD simulations. For T > 4000 K, the aver- 
age band gap that is less than the temperature (dashed line) , 
which implies that the polymeric and the plasma state are 
good electrical conductors. 



of molecular hydrogen. DFT-MD simulations of hydro- 
gen also show that the dissociation of molecules leads to 
a region with dP/dT\v < at high density because hy- 
drogen atoms can also be packed more efficiently than 
molecules [34 | [35 ] . 

When the simulations are cooled from 4000 K, the sys- 
tem remained in a polymeric state that exhibits a lower 
pressure than the original simulations with methane. 
This confirms the free energy calculations in Ref. [55] 
that showed the polymeric state to be thermodynami- 
cally more stable. 

In Fig. |3j we also plot the electronic band gap averaged 
over the trajectory at different temperatures. From 300 
to 3000 K, the system remains in an insulating state with 
a gap of 3 eV or larger. When the system transforms into 
the polymeric state at 4000 K, the gap drops to a value 
close to zero and remains small as the temperature is in- 
creased further and the system transforms into a plasma. 
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One may conclude that this insulator-to-metal transition 
is driven by the disorder that the dissociation introduced 
into the fluid and resembles recent findings for dense he- 
lium P31 [3T] where the motion of the nuclei also led to a 
band gap closure at unexpectedly low pressures. 

If the band gap energy is comparable to fc^T, the sys- 
tem becomes a good electrical conductor and its opti- 
cal reflectivity increases. Since reflectivity measurements 
under shock conditions [32] have now become well es- 
tablished, the transformation into a polymeric state at 
4000 K can directly be detected with experiments. How- 
ever, if this transformation already occurs at lower tem- 
perature, the polymeric transition will be more difficult 
to detect with optical measurements alone because Fig. [3j 
shows that a band gap opens back up as the polymeric 
state is cooled. 
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FIG. 4: Evolution of fraction of single carbon atoms during 
the MD simulation at 4000K and 1.2gcm~ 3 . The decrease 
measures how many methane molecules polymerize to form 
long carbon chains, which is illustrated in Fig. [5] 

We now study the transformation into the polymeric 
state in more detail by analyzing the clusters in the fluid. 
Most simply, one can detect different polymerization re- 
actions by analyzing the carbon atoms only. Using a 
distance cut-off of 1.8 A, we determined how many C n 
clusters were present in each configuration along differ- 
ent MD trajectories. Methane molecules are classified as 
Ci clusters in this approach. Figure [3] shows a drop in 
the fraction of Ci clusters as the system at 4000 K and 
1.2 gcm~ 3 transforms into a polymeric state over the 
course of approximately of 3 ps. 

In Fig. [5j we show a series of snap shots from the same 
DFT-MD simulation in order to illustrate how hydrocar- 
bons chain of increasing length are formed. As soon as 
the first C2 cluster occurs, we also find molecular hydro- 
gen to be present. This resembles typical dehydrogena- 
tion reactions such as nCH 4 — > C n H 2 „+2 + (n — 1)H 2 
that were also analyzed in [H). This implies that molec- 
ular hydrogen may be expelled as a methane rich system 
transforms into a polymeric state. 
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FIG. 5: (color online) Series of snapshots from MD simula- 
tions at 4000 K and 1.2 gem -3 . The large and small spheres 
depict the carbon and hydrogen atoms, respectively. The C- 
C, C-H, and H-H bonds are illustrated by dark thick lines, 
thin lines, and thick light lines, respectively. Panel (a) de- 
picts the initial configuration with 27 CH4 molecules, (b) 
shows the formation of the first ethane molecule. Some H2 
molecules have also formed, (c) and (d) show the first C3 and 
C4 chains that formed, (e-g) illustrate linear C6, C7, and C10 
chains while (h) shows a 12 atoms chain with one branching 
point. 



In the polymeric regime, the system is primarily com- 
posed of different types of hydrocarbon chains but each 
chain exists only for a limited time before it decomposes 
and new chains are formed. Every chemical bond has a 
finite lifetime, which makes it similar to a plasma, but 
many chains are present in every given snap shot, which 
may be a unique feature of carbon rich fluids and makes 
the distinction from a typical plasma state worthwhile. 

We computed the lifetime of each C„ cluster in the 
DFT-MD simulation, neglecting any information that 
may be stored in the position of the hydrogen atoms. 
In Fig. [6] we plot the fraction of carbon atoms that were 
bound in clusters of different size and lifetimes. The plot 
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FIG. 6: Cumulative lifetime distribution of carbon clusters 
of different sizes at various densities and temperatures. The 
square in the upper panel shows that, on average, 26% of the 
carbon atoms are bound in carbon clusters with 20 atoms or 
more that exit for at least 25 fs. Similarly, the open circle 
illustrates that 66% of the atoms are bound in clusters with 5 
atoms or more. Simulations were characterized as polymeric if 
at least 40% of the atoms were bound in clusters with 5 atoms 
or more that existed for at least 20 fs (diamond symbol). 



is cumulative in two respects. The curve for cluster size 
n presents the average fraction of carbon atoms in clus- 
ters with at least n atoms with a lifetime at least as long 
as the lifetime specified on the ordinate. By definition, 
the Ci(t) curve approaches 1 for small t. It is surprising, 
however, that, in simulations at 2.0 gem -3 and 4000 K, 
26% of all carbon atoms were bound in clusters with 20 
atoms or more that existed for at least 25 fs. 

We characterized a simulation as polymeric if, on aver- 
age, at least 40% of the carbon atoms were at any given 
time bound in clusters with 5 atoms or more with life- 
times of at least 20 fs. All simulations with cluster size 5 
curves that fall above the diamond symbol in Fig. [6] are 
then considered to be polymeric. The minimum cluster 
size of 5 atoms distinguishes the polymeric state from 
simulations with stable methane molecules while the life- 
time constraint distinguishes it from a plasma where one 
finds many small ionic species that exist for less than 1 fs. 




Time (fs) 

FIG. 7: Lifetime distribution of plasma species at 20 000 K 
and 1.5gcm -3 . 

Figure [7] shows a typical lifetime distribution of the most 
common ionic species that we obtained from a cluster 
analysis that includes C and H atoms. 
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FIG. 8: Pressure as a function of temperature along the shock 
Hugoniot curves from Fig. [l] for different initial densities, 
po. Shown for reference are the Uranus and Neptune isen- 
tropes [29] . 

This cluster analysis allowed us to characterize all sim- 
ulations in Figs. [T] and [8] as either molecular, polymeric, 
or as a plasma. It is important to note that the clas- 
sification refers to the state reached in our picosecond- 
timescale DFT-MD simulations and, in case where we 
find the CH4 molecules to be stable, this may not neces- 
sarily be the thermodynamic ground state. If all methane 
molecule remained stable in our simulation we described 
the state as molecular, even though our cooled simula- 
tions did not revert back to a molecular state. Previous 
studies (Ref. [32]) showed the polymeric state is thermo- 
dynamically favorable for pressures above 4 GPa. 

In Fig. [5J we plot temperature as a function of pressure 
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along each of the Hugoniot curves. There is a noticeable 
change in the slope at approximately 4000-5000 K as the 
the Hugoniot curves enter the polymeric regime and then 
the plasma state, but it is not as pronounced as indicated 
in the temperature-density plot in Fig. [T] 

We also included in Fig. [8] the isentropes for Uranus 
and Neptune that were approximately determined using 
DFT-MD simulations of only H 2 [55]. These indicate 
that precompressed shock experiments with initial densi- 
ties from 0.35 to 0.70 gem -3 will be capable of reaching 
Uranian or Neptunian interior conditions over a pressure 
range from 20 to 150 GPa. While this pressure range 
can also be accessed with static diamond anvil cell ex- 
periments, the corresponding temperature range of 2000- 
4000 K may be difficult to reach. 

More importantly, Figure [8] shows that the predicted 
planetary isentropes pass through the polymeric regime. 
This implies that any methane ice that was incorporated 
into the interiors of Uranus and Neptune will not have 
remained there in molecular form. Based on our simula- 
tions, we predict it to transform into a polymeric state 
permitting a substantial amount of molecular hydrogen 
to be expelled. The hydrogen gas may be released into 
the outer layer that is rich in molecular hydrogen. 

It is also possible that more complex chemical reactions 
may take place with a dense mixture of water, ammonia, 
and methance ice in the interiors of Uranus and Neptune. 
It could also be difficult to determine the state of chem- 
ical and thermodynamic equilibrium of the ice mixture 
with a single DFT-MD simulation for a given stochiom- 
etry because even very long simulations may remain in 
the metastable state as we have seen here for methane. 

The fact that hydrogen gas may be expelled from hot, 
dense hydrocarbons may also imply that no superionic 
state of methane, in which the carbon atoms remain in 
their lattice positions while the hydrogen atoms move 
through the crystal like a fluid, exists. Superionic be- 
havior has been predicted theoretically for hot, dense 
water and ammonia |33j but no predictions have been 
made for methane despite the fact that the ratio of ionic 
radii of both species are comparable to those in H2O and 
NH3. Since a superionic state requires high pressure and 
temperatures, methane may never transform into such 
a state because hydrogen gas may be released instead. 
The remaining structure may be too dense for hydrogen 
diffusion to occur. 

The structure of the fluid was further investigated by 
examination of pair correlation functions as shown in 
Fig. [9] All pair correlation functions were computed at 
the density 2.0 gem -3 . At 2000 K, methane molecules 
clearly dominate the system. The large peak at 1.05 
A corresponds to the C-H bond of methane. This peak is 
followed by a minimum close to zero, which is indicative 
of stable molecules. When the temperature is raised to 
3000 K, a new C-C peak begins to emerge at 1.4 A, cor- 
responding to the carbon bonds in longer hydrocarbon 
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FIG. 9: Pair correlation functions, g(r), between different nu- 
clei as a function of temperature at a density of 2.0 gem" 3 . 
The intramolecular C-H and H-H peaks at 1.05 and 1.7 A dis- 
appear as CH4 molecules dissociate with increasing tempera- 
ture, while the intensity of the C-C peak at 1.4 A rises with 
increased polymerization and then drops again as a plasma 
state is reached at 10 000 K. 



chains. A corresponding decrease in the C-H peak, and 
the appearance of a H-H peak of molecular hydrogen, is 
also seen. At 4000 K, these processes become more ap- 
parent, as the C-C peak at 1.4 A intensified and the C-H 
peak diminished. The flattening of the H-H curve is a 
consequence of the dissociation of the methane molecules. 
Lastly, at 10 000 K the C-C peak, while reduced, is still 
evident. However, at this temperature all molecules are 
very short lived and unstable, and the broadening of C-H 
and H-H peaks illustrate this fact. 
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FIG. 10: Subset of Hugoniot curves from Fig. [I] for three ini- 
tial densities, po, specified in the legend. The dashed lines 
correspond to polymeric samples that have been cooled from 
simulations at 4000 K. The dash-dotted line refers to simula- 
tions without contributions from excited electronic states. 
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We predicted the shock Hugoniot curves for the poly- 
meric state by cooling the simulations from 4000 to 300 K 
at six different densities from 0.8 to 2.0 gem -3 . Starting 
from the geometries and velocities at 4000 K, we lowered 
the temperature in intervals of 1000 K by rescaling the 
velocities and performing MD simulations for 2 ps at each 
step. The structure of hydrocarbons, that were initially 
present, change very little during the cooling simulations 
because the temperature dropped below the activation 
barrier for such reactions. As expected, all simulations 
eventually froze. For p = 1.2 gem -3 , it attained a com- 
position of C9H20 + C 6 Hi4 + 2C 2 H 6 + 8CH4 + 15H 2 . 
These finding resemble results of the shock experiment 
by Hirai et al. [15] . who observed that the hydrocarbons 
produced did not revert back to methane when the sam- 
ples were cooled. 



Figure 10 shows the Hugoniot curves for a cooled sim- 
ulation corresponding to three different initial densities. 
At 4000 K, they converge to our original Hugoniot curves 
but deviations become increasingly pronounced with de- 
creasing temperature, culminating in a density increase 
of, respectively, 0.23 (31%), 0.21 (24%), and 0.17 gem -3 
(24%) for the initial densities of po = 0.4, 0.5, and 
0.6 gem -3 at 300 K. Because we kept the initial param- 
eters, E and Vq, of methane in Eq. [T] the deviations 
between the original and the polymeric Hugoniot curves 
persist in the limit of zero temperature. We predict the 
findings of a future shock experiment on methane to fall 
in between both sets of curves. For low shock veloci- 
ties, one would expect the results to track our original 
Hugoniot curves for stable molecules. For higher shock 
speeds, the sample in the experiment may convert to a 
polymeric state at lower temperatures than we see in our 
simulations. In this case, we predict a significant increase 
in density to occur that shifts the measured Hugoniot 
curve from the original towards the polymeric Hugoniot 
curve. Figure [3] suggests that this may be accompanied 
by a drop in pressure and a modest change in optical 
properties. 

In shock experiments, one may observe the transfor- 
mation into a polymeric state at lower temperatures than 
we predict with our simulations because experiments last 
for nanoseconds and samples are much larger. However, 
also experiments may also not reach chemical equilibrium 
at lower temperatures because the activitation energies 
for different chemical reactions may be higher than the 
available kinetic energy that could trigger such reaction. 

The fact that the molecular and the polymeric Hugo- 
niot curve diverge below 4000 K in Fig. [10] also has im- 
plications for decaying shock experiments [36] . This new 
experimental technique allows one to map out a whole 
segment of the Hugoniot curve with a single shock wave 
experiment. At the beginning of the shock propagation, 
the material in compressed to a state of high pressure 
and temperature on the Hugoniot curve. As the shock 
wave continues to propagate, the particle velocity is grad- 



ually reduced and new material is compressed to a lower 
pressure and temperature. The shock velocity is contin- 
uously monitored with an interferometer so that a whole 
segment of Hugoniot curve can be measured in one shock 
experiment. If a decaying shock acrosses the polymeric 
regime of methane, we predict the following to occur. As 
the shock enters the polymeric regime, new material at 
the shock front would still be compressed as in a usual 
shock experiments and may yield a state on the molecular 
or on the polymeric Hugoniot, as we discussed. However, 
as material behind the shock front expands from a high- 
temperature state, it will change from a plasma into the 
polymeric state rather than converting back to molecular 
methane, as we have observed in our cooled simulations 
in Fig. [3] This may lead to a structured shock wave with 
various parts in different thermodynamic states. So the 
polymerization of methane may provide a mechanism for 
introducing a shock velocity reversal that were observed 
in [10] and discussed in [TT] . 

Finally, we investigated the effects that thermal elec- 
tronic excitation have on the shape of the Hugoniot 
curves. While electronic excitations were included in all 
simulations discussed so far, we performed some DFT- 
MD simulations where the electrons were kept in the 
ground state. This resulted into drastic reductions of the 
computed pressures and internal energies but both effects 
nearly cancelled [TS] when the Hugoniot curve was cal- 
culated. So Figure [T0| shows only a modest shift towards 
lower densities for temperatures over 10 000 K. 



CONCLUSIONS 

We identified and characterized molecular, polymeric, 
and plasma states in the phase diagram of methane by 
analyzing the of the liquid in DFT-MD simulations that 
we performed in a temperature-density range from 300 
to 75 000 K and 0.8 to 2.5 gem -3 . We presented shock 
Hugoniot curves for initial densities from 0.35 to 0.70 
gem -3 . We predict a drastic increase in density along 
the Hugoniot curves as the sample transforms into a poly- 
meric state. This transformation is accompanied by an 
increase in reflectivity and electrical conductivity because 
the polymeric state is metallic. 

This state is also prone to dehydrogenation reac- 
tions, which has implications for the interiors of Uranus 
and Neptune. These planets' isentropes intersect the 
temperature-pressure conditions of the polymeric regime 
and thus molecular hydrogen may be released from ice 
layers of both planets. This could lead to more com- 
pact cores in Uranus and Neptune, in contrast to recent 
predictions of core erosion for Jupiter and Saturn [2] [3] . 

We also argued that the dehydrogenation reactions 
prevent methane from assuming a superionic state at 
high pressure and temperature that was predicted for 
water and ammonia 133 . However more theoretical and 



experimental work on mixture of planetary ices will be 
required to provide much-needed constraints for the in- 
teriors of ice giant planets [37] . 



[1] W. B. Hubbard, Science 214, 145 (1981). 

[2] H. F. Wilson and B. Militzer. Astrophys. J, 745:54, 2012. 

[3] H. F. Wilson and B. Militzer. Phys. Rev. Lett, 

108:111101, 2012. 
[4] W. J. Borucki et al, The Astrophysical Journal, 726 19 

(2011). 

[5] D. G. Hicks, T. R. Boehly, J. H. Eggert, J. E. Miller, 
P. M. Celliers, and G. W. Collins. Phys. Rev. Lett, 
97:025502, 2006. 

[6] M. D. Knudson and M. P. Desjarlais. Phys. Rev. Lett., 
103:225501, 2009. 

[7] K. K. M. Lee et al. J. Chem Phys., 125:014701, 2006. 

[8] B. Militzer and W. B. Hubbard. AIP Conf. Proc, 
955:1395, 2007. 

[9] R. Jeanloz, P. M. Celliers, G. W. Collins, J. H. Eg- 
gert, K. K. M. Lee, R. S. McWilliams, S. Brygoo, and 
P. Loubeyre. Proc. Nat. Acad. Set., 107:9172, 2007. 
[10] D. K. Spaulding et al. Phys. Rev. Lett, 108:065701, 2012. 
[11] B. Militzer, " Ab Initio Investigation of a Possible Liquid- 
Liquid Phase Transition in MgSiOa at Megabar Pres- 
sures" , submitted to the Journal of High Energy Density 
Physics (2012). 

[12] B. Militzer and D. M. Ceperley. Phys. Rev. Lett, 85:1890, 
2000. 

[13] B. Militzer, D. M. Ceperley, J. D. Kress, J. D. John- 
son, L. A. Collins, and S. Mazevet. Phys. Rev. Lett., 
87:275502, 2001. 

[14] B. Militzer. Phys. Rev. B, 79:155105, 2009. 

[15] B. Militzer. Phys. Rev. Lett., 97:175501, 2006. 

[16] B. Militzer, F. Gygi, and G. Galli. Phys. Rev. Lett, 
91:265503, 2003. 

[17] L. R. Benedetti, J. H. Nguyen, W. A. Caldwell, H. Liu, 
M. Kruger and R. Jeanloz, Science 284 100 (1999) 

[18] H. Hirai, K. Konagai, T. Kawamura, Y. Yamamoto and 
T. Yagi, Phys. Earth. Planet. Inter. 174 242 (1999). 

[19] W. J. Nellis, D. C. Hamilton and A. C. Mitchell, J. Chem. 
Phys. 115 1015 (2001). 

[20] H. B. Radousky, A. C. Mitchell and W. J. Nellis, J. 



Chem. Phys. 93, 8235 (1990). 
[21] F. Ancilotto, G. L. Chiarotti, S. Scandolo, E. Tosatti, 

Science 275, 1288 (1997). 
[22] L. Spanu, D. Donadio, D. Hohl, E. Schwegler and G. 

Galli, Proc. Nat. Acad. Sci. 108 6843 (2011). 
[23] G. Gao, A. R. Oganov, Y. Ma, H. Wang, P. Li, Y. Li, T. 

Iitaka, and G. Zou, J. Chem. Phys 133 144508 (2010). 
[24] N. Goldman, E. J. Reed and L. E. Fried, J. Chem. Phys. 

131, 204103 (2009). 
[25] Y. B. Zeldovich and Y. P. Raizer, Physics of Shock Waves 

and High- Temperature Hydrodynamic Phenomena, Aca- 
demic Press, New York (1966). 
[26] G. Kresse and J. Furthmiiller, Phys. Rev. B 54, 11169 

(1996). 

[27] P. E. Blochl, Phys. Rev. B 50, 17953 (1994). 

[28] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. 

Lett. 77, 3865 (1996). 
[29] R. Redmer, T. R. Mattsson, N. Nettelmann, and M. 

French Icarus 211 798 (2011). 
[30] N. D. Mermin. Phys. Rev., 137:A1441, 1965. 
[31] L. Stixrude and R. Jeanloz. Proc. Nat. Ac. Sci., 

105:11071, 2008. 
[32] P. M. Celliers, P. Loubeyre, J. H. Eggert, S. Brygoo, R. S. 

McWilliams, D. G. Hicks, T. R. Boehly, R. Jeanloz, and 

G. W. Collinsl. Phys. Rev. Lett., 104:184503, 2010. 
[33] C. Cavazzoni, G. L. Chiarotti, S. Scandolo, E. Tosatti, 

M. Bernasconi, and M. Parrinello. Nature, 283:44, 1999. 
[34] B. Militzer, W. H. Hubbard, J. Vorberger, I. Tamblyn 

and S. A. Bonev Astrophys. J. Lett, 688:L45, 2008. 
[35] M. A. Morales, C. Pierleoni, E. Schwegler and D. M. 

Ceperley. Proc. Nat. Acad. Set. 107:12799, 2010. 
[36] J. H. Eggert, D. G. Hicks, P. M. Celliers, D. K. Bradley, 

R. S. McWilliams, R. Jeanloz, J. E. Miller, T. R. Boehly, 

and G. W. Collins. Nature Phys., 6:40, 2010. 
[37] R. Helled, J. D. Anderson, M. Podolak and G. Schubert, 

Astrophys. J. 726:15, 2011. 
Acknowledgements: This work was supported by 
NASA and NSF. Computational resources were supplied 
in part by NCCS and TAC. 



The following table is to be published as online 
supplementary information 



T (K) 


P (GPa) 


J — * l ~ V / v^xxzi. J 




n — fiOO crrm - ^ 
fJ — U.UUU g cm 




300 




— 1 7^/1^1 T 1 




p= U.oUU gem 




300 


O.yo4(lo ) 


r\ i 7om a ( 1 n\ 
— U.l i 2U14(1UJ 


1 nnn 

1UUU 


o.Do^4 ) 


— U. IDOoU^O ) 


9000 

ZrUUU 


i n 4.7/''?^ 

-LU.^ 1 yO ) 


1 c »fi7A/''3 s l 


3000 






4000 


14.37(11) 


— 1 ^^zL^^l 




1 nOO iT/->vn — 3 

p — i.uuu gem 




300 


lo. / o(4J 


n 1^^n ^ >nc;/11^ 
— U.luyzUo(ll ) 


1 nnn 

1UUU 


1 7 QQ/'/l ~\ 

i ( .yo^4 j 


— U. IDZ / 4^ZJ 


2000 


ZU . O^i^ J. -L y 




^000 


94 09/' 1 s ! 

Z4:.UZ^ 1U ^ 


— U. l i ±'±OOI v O ) 


4000 


9^ 00f1 1\ 

£0 . UU ^ J- 


yj. ioz _l \ £ ) 




1 Oni iT/^m - 3 

p — i.zui gem 




oUU 


Zo.ouy (iy j 


— U. 100040^ ) 


1 nnn 

1UUU 


oZ.4Z^0 J 


— U. 10(500^Z ) 


onnn 
zuuu 


oO.O / ^ J-Zj 


— U. 14yZU(^4 J 


Qnnn 
oUUU 


a n n^f i i \ 

4U.U0( ( 11 ) 


— U. loyOl^o ) 


a nnn 

4UUU 


/in /i / q ^ 

4U.4^o ) 


n 1 0*57/' /I ^ 
— U. IZo / ^4J 


OUUU 


4U.41 ^ 14 ) 


— U. 1U0 / [Z) 


£000 
OUUU 


AO fix. /I f\\ 
4o.OO^ lO ) 


— u.uyoui^z j 


7000 

I uuu 


4 i .Ul O J 


O 0S7^0/'1 1 ~\ 
— U.Uo i OU^ll ) 


i nnnn 
1UUUU 


oy. / O^y J 


— U.U04ol^lZ ) 


zuuuu 


1 n7 /i 7/'on\ 

1U i .4 i ^ZU j 


n nnfiQnf t 1 A 
U.UUOoU^ll J 


Qnnnn 
oUUUU 


lOo. 44^10 ) 


n nQ007(^ t n^ 
U.UoZZ ( ylv ) 


zinooo 


91 9 


1 RA^fl 9^ 
U. 1Z 1 


75000 


417.0(3) 


491 50f 13^1 




— J. .OOO g 




2000 


^9 SI H s ! 






p= 1.498 gem ^ 




oUU 


oy.40^0 j 


— u. 10 1 zy4^iu ) 


1 000 
1UUU 


04.00^0 ) 


o i c io ^ ^s/' ^ ^^ 

— U. lOUOO^O ) 


onnn 

ZUUU 


( U.y / J-Oj 


n 1 A C\K7( 
— U. 14U0 '1,0) 


^000 
OUUU 




— U. lOUOO^O ) 


4000 
4UUU 


71 OQYl K\ 

( ± .uy ^ io ) 


11 OP.l'W 
— U. 1 IZOI^O ) 


^000 
OUUU 


i z . ( o ^ i y y 


■ — u.uyo4( ( z ) 


finnn 

ouuu 


77.7(4) 


n noo^ 1 /"i a \ 
— U.Uoy0H i 14 ) 


7000 

I uuu 


oq y/oA 
oo. / y£) 


O 0S1 ^(9} 

— U.UolOI^Z ) 


i nnnn 
1UUUU 


1 nn 7fih /i \ 
1UU. /y\14 j 


— U.U0y41(^o ) 


onnnn 

ZUUUU 


loz. / y^zu ) 


U.UlUoO^ll J 


QOOOO 
OUUUU 


99Q C\((\\ 

zzy .u^o ) 


O 0870^1 1 ^ 

U.Uo 1 U^ll ) 


40000 


9Q9 7(^ 7\ 


o i fi^&n s ! 


7^000 
/ OUUU 








p— 1.600 gem ^ 




9000 

zuuu 


oo.oy ^ id ) 


— u. looyo^o ) 


Qnnn 
oUUU 


nn n x.f~i a\ 
yu.yo^i4 ) 


n 1 OfiQ7M 1 \ 
— U. IZOo 1 


a nnn 

4UUU 


o4.uz^iy ) 


— U. lUoOl^Z J 


c^nnn 
OUUU 


07 cr/oA 

o ( .o{Z ) 


— u.uyooo^io j 


10000 


118.6(2) 


-0.05693(9) 


20000 


185.2(4) 


0.0129(4) 


30000 


254.2(4) 


0.08663(10) 


40000 


326.4(3) 


0.16661(18) 


75000 


599.8(5) 


0.4841(2) 



T fK\ 
1 (KJ 




Hi v / i_ 114 j 




s-\ — 1 77 ^\ cr pm 3 
ft — 1.1 ( J g cm 




9finn 
zuuu 


lie; fxAf] Q\ 


1 1fi1Q( A\ 

— u.iouoy^ j 


^ooo 


1 90 

IZU.O^O 1 


nil QQfif 1 1 ~i 
— u.iiyoo^n j 


aooo 
4uuu 


110 OfzlA 
1 lU.U(4j 


101 eYO 1 ! 
— U. lUlo^Z ^ 


^ooo 

OUUU 


Qfif'l Q\ 
llD.OD^lO ) 


0Q09 c ;/' 1 ^ 
— u.uyuzo^ 10 ) 


i nnnn 
1UUUU 


lOl.oO^lo J 


— u.uozzu^y ) 


90000 
ZUUUU 


99V e^Q\ 
ZZ { .O{o) 


n m VQfi /■ 1 o"\ 
U.U1 1 ooyiz ) 


^onoo 

ouuuu 


OU'i. u ^o y 


0Q1 7((\\ 


aoooo 




1 700^9^1 


7^000 
/ OUUU 




U. 4:000^0 ) 




n — 9 010 trr*m~~ 3 

— z.uiu g uii 




9finn 
zuuu 


1 fil 77^1 7\ 
1U1 • i ((1/ ) 


— u.izuoy^u ) 


3000 


104.4(4 j 


n 1 r\a 1 0^1 £t\ 

— U.lUbl2(luJ 


a nnn 

4UUU 


100. 0(Z j 


— U.UyZ (D[ll) 


OUUU 


i A9 fie fi q\ 
lOZ.Oo^lo J 


n neons/ 1 n^ 
— U.UoZUU^lU ) 


10000 


one; ci/q\ 

zvd.\){o ) 


nn/i/111/1 nA 
— U.U4411(1U J 


90000 

zuuuu 


9Q9 C\(A\ 

zyz.u^4 j 


O 09^0(^9^ 
U-UZOU^Z ^ 


QflOOO 

ouuuu 


0oU.4(0 ) 


u.uyooy ^ 10 ) 


zloooo 


/IRQ f\(fi\ 


1 7^^^^ 
U. 1 1 uo^o ^ 


75000 


SI 1 Offi 1 ) 






n — 9 1 9Q o-r>m — ^ 
p — z.izy gtiii 




1 0000 
-LUUUU 




u.uoyuy^y j 




n — 9 9^7 err*™ - 3 




e^nnn 
OUUU 


OOO r 7('i\ 

ZZZ. (yo) 


n ri70i /i / 1 1 A 
— U.U i Zl4(ll ) 


i nnnn 
1UUUU 


Z f Z.O^O J 


n nQ/i qq/i a \ 
— U.Uo4oo( i 14J 


ZUUUU 


ouy .y (4 j 


u.uo4oy( i ioj 


Qnnnn 
oUUUU 


4oo.y \p ) 


n 1 nvoe/ 1 e^ 
U.lU/ Zo^loJ 


AOOOO 




1 RA^^^l 
u. ±o^±oy^ ) 


75000 


954 7ffi s ) 


f) 4Q1 A7(\ f(\ 




~ — o 07^ a rm — 3 
p — £.0 1 U g cm 




10000 


^08 ^fzH 

OUO . O y 


— 09Q4 C 1(' 1 9~l 




p — z.ouz gem 




OUUU 


9Q9 A f 

zyz. 4^o ) 


OKI oe/ 1 Q^ 
— U.UOlUo^loJ 


10000 


349.7(5) 


-0.02343(16) 


20000 


457.5(2) 


0.04569(10) 


30000 


566.9(9) 


0.1180(2) 


40000 


679.6(5) 


0.19499(14) 


75000 


883.4(7) 


0.32519(17) 



TABLE I: Pressure and internal energy from our DFT-MD 
simulations at various densities and temperatures. The 1 a 
error bars are given in brackets. 



